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Critical behavior of solitonic waveforms of Bose-Einstein condensates in optical lattices (OL) 
has been studied in the framework of continuous mean-field equation. In 2D and 3D OLs bright 
matter-wave solitons undergo abrupt delocalization as the strength of the OL is decreased below 
some critical value. Similar delocalizing transition happens when the coefficient of nonlinearity 
crosses the critical value. Contrarily, bright solitons in ID OLs retain their integrity over the whole 
range of parameter variations. The interpretation of the phenomenon in terms of quantum bound 
states in the effective potential is proposed. 
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I. INTRODUCTION 

Bose-Einstein condensates (BEC) display remarkably 
rich nonlinear wave phenomena, among which solitons 
are particularly interesting. To date both common types 
of solitons, bright and dark, are experimentally realized in 
BECs, respectively, with attractive and repulsive inter- 
atomic forces [1,2]. Theoretical studies of the properties 
of matter- wave solitons has resulted in successful descrip- 
tion of basic features of soliton dynamics in BEC [3]. 

Considerable effort is being devoted to investigation 
of BEC in periodic potentials created by laser standing 
waves - so called optical lattices (OL). Subjecting a BEC 
to a periodic potential leads to a number of interesting 
phenomena, including the formation of a specific class 
of localized modes called gap solitons. The definition 
'gap' refers to the fact that this type of solitons exist in 
the band gaps of the matter-wave spectrum. Conceptu- 
ally important point here is that, the gap solitons can be 
composed of repulsive atoms, which makes them favor- 
able against bright solitons in continuous attractive BEC 
tending to collapse at high atomic densities. 

After the first discussion of the possibility of bright 
solitons in repulsive BEC loaded in OL [4], substantial 
theoretical research was performed aimed at understand- 
ing their development and properties [5]. One of the 
recent advances in this direction has been the proof of 
the possibility of matter-wave solitons in two- and three- 
dimensional (2D and 3D) OLs. Particularly, the emer- 
gence of multi-dimensional gap solitons in repulsive BEC 
arrays due to the phenomenon of modulational instabil- 
ity was demonstrated in [6]. Different modes of stable 
2D gap solitons in OLs were shown to exist by numerical 
solution of the corresponding 2D Gross-Pitaevskii equa- 
tion [7]. Multidimensional solitons in media governed 
by the self-focusing cubic nonlinear Schrodinger equation 
(NLSE) with a periodic potential, of which an attractive 
BEC in OL is a particular example, are studied in [8] by 
means of the variational approximation and direct nu- 
merical simulations. 

While have not been experimentally realized yet, 
bright matter-wave solitons in OLs currently are among 



the intensively explored subjects in BEC. An important 
aspect of OLs is that they provide a unique possibility 
to investigate the behavior of localized excitations and 
solitons in both continuous and discrete ends of the spec- 
trum by adjusting the strength of the OL. A fundamental 
issue that the cubic NLSE does not support stable soli- 
tons in dimensions higher than one, while its discrete 
analog does, can be addressed more effectively involving 
the BEC in OL. One relevant problem - the delocalizing 
transition of BECs in deep 2D OLs has recently been in- 
vestigated in the framework of discrete NLSE [9]. The 
analysis of Ref. [9] is relevant to the tight-binding approx- 
imation, when the BECs in neighboring lattice sites are 
weakly linked. This approach also relies on the possibil- 
ity, that the tunneling amplitude between OL sites can be 
accurately determined. In the case of uniformly loaded 
OL the tunneling amplitude can be estimated using the 
Wannier function basis [10]. However, in the case of lo- 
calized excitations involving just few lattice sites with 
significantly different occupation numbers, calculation of 
the tunneling amplitude becomes difficult. 

The aim of this paper is twofold. From one side we 
present a detailed numerical investigation of the delocal- 
izing transition beyond the tight-binding approximation, 
when the system is described by the continuous mean- 
field equations. To this regard we provide both analyt- 
ical and numerical evidences for the existence of stable 
multidimensional solitons in OLs, and we study the fade 
and recovery of the solitonic structures as the strength 
of the optical lattice or coefficient of nonlinearity adia- 
batically varies in time. The major topic of interest will 
be a critical behavior of solitonic waveforms in a periodic 
potential, expressed as sudden disintegration of the soli- 
ton when the system parameters attain particular values. 
This will be done for BECs with both positive and nega- 
tive scattering lengths. From this study it emerges that 
the variational ansatz while correctly predicting the sta- 
bilization of solitons due to the OL, it fails to predict the 
existence of a delocalizing transition. 

From the other side, we address to the physical mech- 
anism underlying the delocalizing transition. To this re- 
gard we propose a quantum interpretation of the phe- 
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nomenon, by showing that the dclocaUzing transition is 
associated to the disappearance of bound states from 
the sohton effective potential in the underlying periodic 
Schrodinger equation. Since in one dimension any po- 
tential well can support at least one bound state (this is 
true oven for infinitesimal well depths), the above inter- 
pretation automatically implies that ID solitons cannot 
undergo delocalizing transitions. This is indeed what we 
find for the ID case. In particular wc show that by de- 
creasing the strength of the OL or the nonlinearity in the 
system, the soliton becomes very extended (in the limit of 
zero nonlinearity it reduces to a Bloch state at the edge of 
the Brillouin zone) , but always recovers its original shape 
when the parameters arc reversed. In contrast, in 2D and 
in 3D cases, we find that there is always a critical value 
of system parameters below which the recovering of the 
soliton becomes impossible i.e. the soliton irreversibly 
disintegrates into extended states. By approximating the 
soliton effective potential with a suitable solvable poten- 
tial and by adopting a variational ansatz description for 
the soliton, we show that the above bound state interpre- 
tation leads to an analytical expression for the occurrence 
of the delocalizing transition which is in good agreement 
with direct numerical simulations of the Gross-Pitaevskii 
equation, this providing support for our approach. 

We remark that the proposed mechanism for the de- 
localizing transition is general and expected to be valid 
also for nmltidimensional intrinsic localized modes or dis- 
crete breathes of nonlinear lattices, as well as for multi- 
dimensional soliton solutions of the NLSE with nonlocal 
interactions. 

Finally, we have explored the role of dimensionality 
by comparison of results for 2D and 3D OLs. The ad- 
vantage of the present approach is that, all the system 
variables arc readily connected to actual physical param- 
eters (e.g. the OL strength to the laser intensity and 
the coefficient of nonlinearity to the atomic scattering 
length, etc.), which makes the experimental verification 
more feasible. 

The paper is organized as follows. In Sec. II we in- 
troduce the model and derive variational equations for 
soliton parameters. Then we check the validity of VA 
equations by comparing with direct numerical solution 
of the Gross-Pitaevskii equation. In Sec. Ill we pro- 
pose a physical mechanism associating the existence of 
solitons and their delocalizing transition with quantum 
bound states in the effective potential. The existence re- 
gion for 2D solitons is presented. In Sec. IV we discuss 
the possibility of experimental observation of the delocal- 
izing transition of solitons. Finally, in Sec. V we conclude 
the results of this study. 

II. THE EXISTENCE OF SOLITONS IN 
PERIODIC POTENTIALS 

The existence, stability and some other properties of 
solitons in ID periodic potentials are explored in the 



context of repulsive BEG in OLs [4,5]. Solitons in self- 
focusing NLSE with a ID periodic potential, with rele- 
vance to particular problems of nonlinear optics is consid- 
ered in [11]. While the subject is well understood in ID 
case, multidimensional solitons in periodic potentials are 
in early stage of research [6-8] . Below we analyze the ex- 
istence of multidimensional solitons in periodic potentials 
in the framework of the variational approximation (VA) 
for both attractive and repulsive interactions. The results 
will then be compared with numerical simulations of the 
full problem. We remark that the VA analysis based on 
a Gaussian waveform was recently considered in Ref. [8] 
for the case of 2D BEG with attractive interactions. In 
this case it was shown that bright solitons are well ap- 
proximated by Gaussian profiles, this being particularly 
true when all the matter is confined into a single cell of 
the OL (single cell solitons). In the case of BEG with 
repulsive interactions, however, satellite peaks appear in 
the soliton waveform and the VA based on a Gaussian 
ansatz leads to wrong results. In this case we will show 
that a Fraunhofer diffraction pattern is a better ansatz 
for the soliton shape, leading to VA equations which are 
in good agreement with numerical simulations. 

A. Basic equation 

The analysis will be based on the mean-field Gross- 
Pitaevskii equation (GPE) for the macroscopic wave 
function of the condensate tp{f,t), which describes the 
properties of a dilute gas near-zero temperature BEG [12] 

in^=[-^V' + Ve,t{f,t)+g\i;\']i,, (1) 

where m is the atomic mass, g is the nonlinear cou- 
pling parameter related to the s-wave scattering length 
as through g = Anh^as/m, Vext{f,t) is the external po- 
tential, generally consisting of harmonic confinement plus 
periodic potential of the OL 

Vext{f, t) = ^ujjrj + s cos{2kLri), (2) 

where = x, y, z and sum over i is assumed. We 
shall be interested in the behavior of localized excita- 
tions, occupying few central unit cells of the OL, when 
the effect of confining parabolic potential is not essen- 
tial. Therefore, in what follows we assume the presence 
of the OL only. For convenience we rescale the Eq.(l) 
by introducing the dimensionless variables — > k^ri, 
t {hkl)/{2m) ■t,e^ e/Er, ip ip/y^- Where 
fcz, = 27r/A, Er = {lT?k\)/{2m), X,kL being the laser 
wavelength and wave vector correspondingly, no is the 
density of BEG. We also designate the coefficient of non- 
linearity as X = 87rnoOs/fc£. Then we have Eq.(l) in the 
form 

i^ = _^2^ + Voi{r,t)i, + x\i'f^, (3) 
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where Voi{r,t) = e[cos(2a;) + cos(2j/) + cos(22;)] is the OL 
potential. 

B. Variational approach 

To deal with compact expressions, we explicitly con- 
sider the 2D case, since the extension to 3D is straight- 
forward. We look for the stationary solution of Eq.(3) in 
the form 

il){x, y,t) = U {x, y) exp{-ifit), (4) 

where /i is the chemical potential. Then the follow- 
ing equation determines the time-independent solution 
U{x,y) 

\/^U + fiU + Voi{x,y}U + xU^ = 0. (5) 

Following the standard procedure of the VA [13,14], we 
consider the effective Lagrangian 

1 f°° 1 
L=^J [mf - i^U' - Voi{x,y)U^ - -xU^]dxdy 

(6) 

which generates the Eq.(5) by minimization 5L = 0. The 
next step of the variational approach involves the choice 
of a suitable trial function - ansatz for the solution. As 
we shall see later, solitonic waveforms in attractive and 
repulsive BECs in the periodic potential have different 
spatial features. This implies the choice of different trial 
functions, which we shall consider separately. 

1. Attractive BEC 

The nature of interatomic forces in BEC can be at- 
tractive or repulsive, which leads to very different prop- 
erties of corresponding BECs [12]. The attractive case 
corresponds to the negative sign of the s-wave scattering 
length as, and the underlying GPE has a self- focusing 
nonlinearity (x < 0). The collapse of BEC at some criti- 
cal number of atoms N ~ 1400 is the main consequence 
of the self- focusing nonlinearity [12]. The presence of a 
periodic potential, however, significantly changes the sit- 
uation and leads to stable localized states. In moderately 
strong and weak periodic potentials (which is the case at 
delocalizing transition point considered below) multidi- 
mensional solitons have a single cell structure and are 
well described by the VA with Gaussian ansatz [8] 

C/(x,y) = ^exp[-|(x^+y2)], (7) 

where the amplitude A and width parameter a are fixed 
by the total number of atoms in thc> condensate = 
\U{x,y)\^dxdy = ttA^ /a. Variational equations are 
derived by substituting the ansatz (7) into the effective 



Lagrangian (6), performing the spatial integration and 
subsequent minimization with respect to variables A and 
a. The equations associate the parameters of the soliton 
with the total number of atoms, strength of the OL, and 
chemical potential of the BEC 

Ar=!!Ifi_^e-V«), ^ = ^(2-a)e-V«-a. (8) 
X V a / a 

A noteworthy property of Eq.(8) is that, for a given 
strength of the OL e there exists a minimal norm Nmin = 
^(1 — 8ee~^) attained at a = 1/2, below which the lo- 
calized solutions do not exist, as illustrated in Fig.l. The 
threshold vanishes for stronger OLs, e > Scr = e^/8 = 
0.92. The existence of solitons in 2D attractive BEC from 
above is limited by the onset of collapse at the critical 
norm Ncoi = 47r predicted by VA [15] (the exact value 
being Ncoi = 11.6993 [16]). 




a 

FIG. 1. The Eqs.(8) predict two localized solutions for a 
given norm A^, one of which (o > 1/2) is stable and the other 
(o < 1/2) is unstable. 

The VA predicts two localized states with different 
widths for a given norm Nmin < N < Ncoi- The stability 
of these two solutions can be inferred from the Vakhitov- 
Kolokolov criterion [17]. According to this criterion, the 
branch of solutions with a > 1/2 has the negative slope 
dN/dn < 0, and therefore is stable. The other branch 
(a < 1/2) appears to be unstable. 



2. Repulsive BEC 

The main feature of a soliton of repulsive BEC (x > 0) 
in a periodic potential is that, it is composed of a cen- 
tral peak and symmetrically spaced satellites residing in 
neighboring cells of the OL. The issue of vital importance 
for VA is the adoption of a suitable ansatz. A Gaussian 
waveform obviously is not satisfactory in this case since it 
doesn't take into account the composite shape of solitons. 
We introduce the following ansatz 
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U{x,y)=A- 



sin (ax) sin (ay) 



ax 



ay 



(9) 



which accounts for the sateUites, and has the advantage 
of yielding simple VA equations. Analytical transforma- 
tions similar to the previous case leads to following set 
of equations 



N 



The relation connecting the norm and soliton parameters 
in this case is N = ■k'^ / a^ , which yields in combination 
with Eqs.(lO) 



A 



97V3/2£ 



27r(iVx + 37r2) 



1/3 



(11) 



Unlike the attractive case, the ansatz Eci.(9) yields a sin- 
gle set of soliton parameters a, A for a given norm N and 
strength of the OL e. In the next subsection we compare 
the prediction of the VA with direct numerical solution 
of the GPE (3). 
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FIG. 2. The Gaussian waveform Eq.(7) with parameters 
prescribed by VA A = l.Q, a = 1.3, N = 2n, inserted as ini- 
tial condition to Eq.(3) evolves into a stable 2D soliton, per- 
forming transient oscillations. The strength of OL is e = 0.92. 
a) Comparison VA vs. PDE, b) stable 2D soliton of attractive 
BEG. 



C. Comparison of VA and PDE simulations 

Validity of the VA is usually verified by comparison 
with the results of direct numerical solution of the un- 
derlying PDE. Below we solve the Eq.(3) by split-step 
procedure using the multi-dimensional fast FoTiricr trans- 
form [18] on the grids of 512, 256x256 and 128x128x128, 
respectively in ID, 2D and 3D cases. The domain size 
and time step are x,y,z £ [— 47r,47r], 5t = 0.001. The 
essential point is the use of absorption on the domain 
boundaries to prevent re-entering of linear waves emitted 
by the soliton during its formation (or under perturba- 
tions) into the integration area. 

The comparison between VA and PDE simulations can 
be done in different ways. We consider the situations 
more relevant to experiments, namely, we change the 
strength of the OL e or the coefficient of nonlincarity x 
(which is equivalent to changing the norm at constant 
x), and follow the evolution of the soliton's amplitude 
according to PDE (3). Then we correlate the numeri- 
cally obtained dependence A{e) or A{N) with analytical 
prediction of VA Eqs.(8,10). Below we perform the com- 
parisons for 2D case. 

1. Attractive case 

At first we need to generate a stable 2D soliton of 
Eq.(3). For this we insert the Gaussian waveform (7) 
with parameters specified by VA (see Fig.l) into the 
Eq.(3) as the initial condition. The waveform undergoes 
some evolution and attains the stable state as shown in 
Fig.2. 



The stability of multidimensional solitons, being a 
very important issue, is worth of special investigation. 
Here we refer to our numerical tests, which show that a 
small deformation of the soliton causes weak oscillations 
around the stable state. After sufficiently long time the 
waveform relaxes back to its original shape proving that 
the stable state is a fixed point, i.e. the soliton is linearly 
stable. The tests also involve the long-time propagation 
of a fundamental soliton with PDE for a period much 
exceeding the characteristic time of the problem. 



After the stable 2D soliton is created, we use it as 
initial condition for Eq.(3) and follow the further evolu- 
tion under time-dependent parameters s{t) or x(0- 
important point to be stressed here is the condition of 
adiabaticity in variation of parameters. In the present 
context, the adiabaticity refers to the situation, when the 
variation of parameters does not excite collective modes 
of the condensate, which means that the ramp times 
should satisfy the condition tramp > fi/lJ-- In our dimen- 
sionless units the above condition implies tramp > 50 for 
a typical condensate with chemical potential of /z = 200 
Hz. The soliton's amplitude as a function of the coeffi- 
cient of nonlincarity, obtained from PDE simulations and 
compared with the prediction of VA is presented in Fig. 3. 
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FIG. 3. The amplitude of a 2D soliton of attractive BEG 
as a function of the coefficient of nonlinearity, obtained by 
numerical solution of the GPE (1) (solid line), and as pre- 
dicted by VA Eq.(8) (dashed line). VA adequately describes 
the dependence ^(x) above the delocalizing transition point 
at X = 0.38. 

An important remark concerning the Fig. 3 is that 
the VA fails to account for the delocalization of a soliton 
occurring when the norm drops below some critical value 
(x = 0.38, for e = 0.92), which is manifested as a rapid 
spreading of the waveform. 



2. Repulsive case 

Similarly to the previous case, wc generate a stable 2D 
soliton of repulsive BEC inserting the waveform Eq.(9) 
with A = 2.5, a = 1.1 as initial condition for the Eq.(3) 
with e = 4.0 and x = 1- After some transition period a 
stable soliton is formed as shown in Fig. 4. In this case 
the oscillations of the amplitude around the stable state 
is less pronounced due to repulsive nature of the conden- 
sate. 




FIG. 4. Formation of a 2D soliton from the initial state 
given by Eq.(9) with A = 2.5, a = 1.1 in the OL of strength 
e — 4.0. a) Transient oscillations of the amplitude, b) stable 
2D soliton of repulsive BEC. The norm has decreased by 50 
% during the formation of a stable soliton. 



The fact that the ansatz (9) properly accounts for the 
satellites can be seen from the comparison of contour 
plots in Fig. 5. 
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FIG. 5. Contour plots of the ansatz Eq.(9) (a), and the 
stable 2D soliton (b) of repulsive BEC. 

Fig. 6 illustrates the soliton's amplitude function 
of the strength of the OL given by VA Eq.(ll), and ob- 
tained by direct numerical solution of Eq.(3). 
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FIG. 6. The amplitude of a 2D soliton of repulsive BEC as 
a function of the strength of the periodic potential. Solid line 
- numerical solution of Eq.(3), dashed line - prediction of VA 
Eq.(ll). 

Here again, as in the attractive case, we see that the VA 
fails to predict the delocalizing transition (at e ~ 1.2). 
Below we consider in detail the delocalizing transition of 
multidimensional solitons in periodic potentials. 



III. DELOCALIZING TRANSITION 

The delocalizing transition of solitons is manifested as 
irreversible transformation of the localized waveform to 
the extended state. The transition can be induced by 
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decreasing of the strength of the periodic potential, or 
reducing the coefficient of nonhnearity in Eq.(3), below 
some critical value. As pointed out in the previous sec- 
tion, the delocalizing transition of solitons is missed in 
VA description. In this section we present an interpre- 
tation of the phenomenon, involving the quantum bound 
states in the effective potential created by the soliton. 
For this we assume the condensate to have the form of 
a single cell soliton, and consider the VA equations valid 
until the delocalizing transition occurs. An example of 
this is shown in Fig. 3 for the case of attractive interac- 
tions. The underlying idea is to consider the GPE (3) as 
a linear Schrodinger equation with the effective potential 

Veff^Vol+x\M^ (12) 

where Vqi is the periodic potential of the OL, and 
the second term represents the contribution of the soli- 
ton. Finding the stationary solution of the Schrodinger 
equation with the potential (12) is similar to the self- 
consistent Hartree-Fock approximation of quantum me- 
chanics. From this point of view it is clear that the exis- 
tence of matter- wave solitons is linked to the existence of 
quantum bound states in the effective potential (12). We 
remark that the periodic part Vqi can be eliminated from 
the effective potential. This is particularly connected to 
the possibility, that the GPE (1) with a periodic poten- 
tial can be reduced to the standard NLSE in the effective 
mass formalism [19]. In other words, the contribution of 
Vol in the linear Schrodinger problem can be accounted 
in terms of an effective mass description for the extended 
states of the linear problem. These states have positive or 
negative effective masses depending on their interaction 
being, respectively, repulsive or attractive (they resem- 
ble electrons or holes of usual solids). In both cases we 
reveal that the effective potential reduces only to the soli- 
ton part with renormalized coefficient of nonlinearity x 
[19] 

Veff = -\xmx,y)f, (13) 

where tpi^x, y) denotes the solution of the nonlinear prob- 
lem in the presence of OL, evaluated at the delocalizing 
transition point. Notice that the effective potential is al- 
ways negative, i.e. the soliton acts always as a potential 
well in the corresponding Schrodinger problem (in the re- 
pulsive case the sign of the potential is reversed by the 
negative effective mass). The problem then is to investi- 
gate the existence of states with negative energy of the 
following Schrodinger equation 

{-V^ + V,fi-E)ct>{x,y) = Q. (14) 

For this one can resort to numerical schemes for Eq. (14) 
to find that at the delocalization point the effective po- 
tential has the critical strength to support just a single 
bound state. In the following, however, we shall pro- 
vide an analytical estimation of the number of bound 
states employing the VA solutions discussed in the pre- 
vious sections and approximating the effective potential 



with a solvable potential for which this number is exactly 
known. To this end we use the following Poschl- Teller 
potential [20] 

VpT{r) = (15) 

cosh(rv2oln2)^ 

to approximate Vef f at the delocalization transition. 
The parametric form of Vpt is fixed so that Vpt 
has the same amplitude and the same integral value: 
2ti rVefj{r)dr, as for Vejj- In Fig. 7 Vef j and its 
Poschl- Teller approximation are depicted for the param- 
eter values corresponding to the delocalization transition 
of Fig.3. 
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FIG. 7. The effective potential of a 2D soliton of attrac- 
tive BEG at the delocalizing transition point, obtained from 
numerical solution of the GPE (3). The dashed line refers to 
the Poschl- Teller potential in Eq. (15), while the dash-dotted 
line is the Gaussian approximation Eq.(7) with parameters 
Ad = 1.38, Od = 1.25, corresponding to the delocalizing tran- 
sition point of Fig.3 at x = 0.38 for e = 0.92. 

The number of bound states of the Poschl- Teller po- 
tential is given by the integer part of n [21] 

This equation can be reduced to a more convenient form 
by using the expressions for the amplitude A and the 
width a of \\iq soliton as obtainc^d from VA. At the de- 
localizing transition n = 1 so that, after substitution 
A^/a = N/tt, the Eq.(16) reduces to ^ = In 2. This 
means that there exists a critical value Ncr = 47rln2/x 
for the number of atoms, below which the soliton dis- 
integrates (notice that this delocalization threshold is a 
factor of In 2 smaller than that of unstable 2D Townes 
soliton [16]). Since Ncr is related to the amplitude and 
width Ad, fld, of the soliton at the delocalization point by 
the variational equations, we can express the delocaliza- 
tion condition in terms of e, a^, A^, as 

e=ia^(l-ln2)ei/-^ ^=1^2. (17) 

In the following we shall compare these expressions for 
the delocalization transition with the results of direct 
PDE simulations. 
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A. ID optical lattice 



Since the ID Schrodinger equation has bound state so- 
lution in any confining potential (even of infinitesimal 
depth), this implies that the delocalizing transition of 
bright solitons cannot occur in ID case. It is interest- 
ing to start with consideration of this situation, since it 
provides more insight into the general problem. 

The dynamics of the condensate is governed by ID ver- 
sion of the GPE (3). The stationary localized solution to 
this equation can be found by standard numerical relax- 
ation procedure. One of such gap soliton solutions for 
repulsive case (x > 0) is presented in Fig. 8 




- 



FIG. 9. When the strength of the OL eo = 1 is decreased 
to zero at t = 50 and then increased back to its initial value, 
the ID gap soliton of repulsive BEG fully restores its original 
form. 




FIG. 8. The initial waveform corresponding to stationary 
solution of Eq.(3) with parameters e = 1, x = 1- Dashed line 
represents ID OL potential V{x) = ecos(2a;). 



Now we explore the behavior of ID soliton under adi- 
abatic variation of the strength of OL e{t) or coefficient 
of nonlinearity xC^); using the above solution (Fig. 8) as 
initial condition for time-dependent GPE (3). We as- 
sume the simplest linear function for variation of the 
amplitude s{t) = so ■ f{t) and coeflacient of nonlinear- 
ity X{t) = Xo ■ f{t), with 



fit) 



2a{t/tend), 



if < i < 0.5 • tend 
if 0.5 • tend <t < tend 

(18) 



The opposite situation also reveals the solitarily of the 
waveform (Fig. 8). Specifically, when the strength of the 
OL is increased enough the BEG matter is entirely pulled 
into the central unit cell (Figs. 10, II). 




50 



FIG. 10. All the BEG matter is pulled into the central unit 
cell as the strength of the OL is linearly increased from e = 1 
to e = 5 during t = 50. 

The following Fig. 11 represents the waveform in OL 
potential at t = 50. 



In Fig. 9 we report the transformation of the solitonic 
state into the extended one, and restoring its original 

shape as the strength of the OL is adiabatically decreased 
to zero s.tt = 50, and then increased back to initial value 
at t = 100. 



FIG. 11. A ID gap soliton of Fig. 8 contracts into the cen- 
tral unit cell when the strength of the OL is increased from 
£ = 1 to e = 5. Dashed line represents cos(2x) 

In the case of attractive BEG (x < 0), the stationary 
solution of the GPE (3) can be found similarly to above, 
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by relaxation method. Below we consider the solution 
corresponding to the strength of the OL e = 1.0, and the 
soliton initially fits a single cell of the periodic potential. 
Then we use this solution as initial condition for GPE (3) 
and observe the spreading and contracting of the wave- 
form as the coefficient of nonlinearity is decreased to zero 
and increased back to its original value, as displayed in 
Fig.l2 




FIG. 12. Soliton of attractive BEG in ID OL of strength 
£ = 1 retains its integrity while the coefficient of nonlinearity 
is decreased to zero at t = 500, and then increased back to 
the original value x = — 1 at t = 1000. 

To conclude this subsection we note, that the bright 
solitons in ID periodic potentials do not exhibit the de- 
localizing transition as the OL is weakened or the co- 
efficient of nonlinearity is decreased. The integrity of 
solitons is manifested by recovering their original shape 
as the periodic potential or coefficient of nonlinearity is 
restored to the initial value. This is in agreement with 
the proposed physical mechanism associating the delo- 
calizing transition with quantum bound states in the ef- 
fective potential created by the soliton. As noted above, 
in ID there are always bound states in the confining po- 
tential, and that is why the soliton doesn't disintegrate 
even when the depth of the effective potential becomes 
infinitesimal. 



B. 2D optical lattice 

We separately consider the delocalizing transition of 
bright solitons of attractive and repulsive BECs in 2D 
OLs. Starting point in thisc cases, similarly to the pre- 
vious situation, is construction of a stationary solitonic 
state in a 2D OL. 



1. Attractive case 

The stationary soliton of attractive BEG in 2D OL 
(Fig. 2), generated from the Gaussian waveform with pa- 



rameters prescribed by the VA was employed as initial 
condition in the GPE (3) with time-dependent param- 
eters s{t), or x(f). By adiabatic variation of these pa- 
rameters one can bring the soliton close to the delocal- 
izing transition point, and then return back recovering 
the initial waveform, or induce the delocalizing transi- 
tion by slightly more decreasing the critical parameter, 
as demonstrated in Fig. 13 by numerical solution of the 
GPE (3). 




FIG. 13. A 2D soliton of attractive BEG recovers its origi- 
nal amplitude when the parameter x{t) is decreased and then 

increased back without crossing the critical value (near point 
G), and irreversibly disintegrates when crossed (point D). A 
small deficit in the return value of the amplitude is caused by 
the energy loss due to imperfect adiabaticity of the process. 



Repeating similar delocalizing transition simulations 
one can establish the existence region of 2D solitons in 
the parameter space e vs. x- I^i Fig. 14 we present the re- 
sult of such numerical simulations for different strengths 
of the OL. Unlike the Fig. 13, the coefficient of nonlin- 
earity x{t) was linearly decreased until zero, as we are 
interested in the critical values Xd leading to delocaliza- 
tion of the soliton. 
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FIG. 14. Delocalizing transition of a 2D soliton in OLs 
of different strength e (shown to the right of corresponding 
curves) . 
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Fig. 15 represents the existence region for 2D solitons 

of attractive BEG in tlie OL. 




0.0 I ■ ' ■ ' ■ ' ■ 1 

0.5 0.6 0.7 0.8 0.9 

X 

FIG. 15. The existence region of 2D solitons of attractive 
BEG. Solid line - parametric solution of Eq.(17), also using the 
conserved norm A'' = TvA^/ad- Stars - delocaliziug trausition 
points of Fig. 14. The inset shows the number of quantum 
bound states evaluated at the delocalizing transition points 
according to Eq.(16). 

The number of quantum bound states in the effective 
potential at the point of delocalizing transition, evalu- 
ated from numerical simulations (Fig. 14) and using the 
Eq.(16), appears to be very close to one (see the inset in 
Fig. 15) for all values of s, which supports the proposed 
physical model. 



2. Repulsive case 



As pointed out, matter- wave solitons of repulsive BEG 
in OLs have a composite structure. The satellites are 
more pronounced when the soliton is driven close to the 
delocalizing transition point in the parametric space s vs. 
X, as illustrated in Fig. 16. 





FIG. 16. a) A stable 2D soliton of repulsive BEG (x = 1) 
in OL of strength e = 4.5. b) Tiic satellites arc more pro- 
nounced when the delocalizing transition point is approached 
by reducing the coefficient of nonlinearity until x = 0.4. 

Similarly to attractive case, the soliton can be brought 
close to the delocalizing transition point by adiabatic 
change of parameters e or %. Then, if the initial values of 
parameters are restored without crossing the threshold, 
the soliton recovers its original shape, otherwise it irre- 
versibly disintegrates. Such a numerical simulation with 
a soliton of repulsive BEG is presented in Fig. 17. 

..Ji (a) 



0.6 




FIG. 17. Diagonal section profile for the time evolution of 
a 2D soliton according to GPE (3). (a) The waveform re- 
covers its original shape when the magnitude of the periodic 
potential is linearly decreased from eo = 4.5 until Smin = 3.42 
at t = 500, and then increased back to co at t = 1000. (b) 
Abrupt delocalizing transition occurs as the strength of the 
OL is lowered below the critical value Sc = 3.38 at t = 500. 
The Eq.(18) for variation of e{t) = eo ■ fit) is applied with 
= 1000, and a = 0.25 for (a), a = 0.26 for (b). 

Of particular interest is the interplay between the cen- 
tral peak and satellites of the lattice soliton in repulsive 
BEG. Some information can be obtained by recording the 



9 



amount of BEC matter I{t) in the central peak (confined 
to a unit cell), while the parameters £ or x adiabatically 
varied in time. 



lit) 



7r/2 



7r/2 



\i^{x,y,t)\'^dxdy. 



(19) 



Time evolution of this quantity is presented in Fig. 18 
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FIG. 18. Amount of BEC matter in the central peak of a 
2D soliton according to Eq.(19), corresponding to Fig. 17b. 



From this figure one can judge about the role of the 
central peak in the integrity of the soliton. The delocal- 
izing transition of a composite soliton implies the rapid 
spreading of its central peak. 



FIG. 19. Wave profile along the main diagonal of a cubic 
domain x,y,z G [— 47r,47r] as obtained by numerical solution 
of the GPE (3) with e{t) = eo • (1 — t/tend) (a), and corre- 
sponding integral Eq.(19) over the central unit cell (b), for a 
3D OL with £0 = 4.5, x = 1-0, tend = 1000. 

Also note the absence of transient oscillations of the 
integral atomic number after delocalizing transition in 
Fig. 19b, as opposed to 2D case (Fig. 18). This is another 
effect of higher dimensionality. 



C. 3D optical lattice 



There is a growing interest in properties of BEC in 
3D OLs [22,23]. The effect of dimensionality in the pro- 
cess of delocalizing transition of matter-wave solitons is 
of particular interest. Since we consider the phenomenon 
as rapid spreading of the wave-packet through tunneling 
of BEC into neighboring lattice sites, the main mani- 
festation of the dimensionality should be the increased 
sharpness of the transition (because of more neighboring 
cells compared to 2D case). 

To verify the above conjecture, we performed numeri- 
cal simulation of the delocalizing transition of a soliton in 
3D OL. The procedure for preparation of a stable soliton 
and subsequent adiabatic variation of system parameters 
are similar to 2D case. The result is presented in Fig. 19. 
As expected, the delocalizing transition in 3D appears to 
be more sharp compared to 2D case. 



IV. EXPERIMENTAL FEASIBILITY 

Addressing the possibility of experimental observation 
of phenomena considered in this paper, we note the rapid 
progress in manipulation techniques for BEC in OLs 
[22-24]. However, the matter- wave solitons in periodic 
potentials of OLs have not been experimentally realized 
yet. Different approaches were proposed, of which par- 
ticular interest is the generation of gap solitons in repul- 
sive BEC employing the phenomenon of modulational 
instability [5,6]. Another possibility is the independent 
successive formation of BEC in single sites of the OL as 
proposed in Ref. [26]. The BEC initially occupying a sin- 
gle or few lattice sites of the periodic potential quickly 
transforms into the solitonic waveform (certainly, if the 
parameters are in the existence region) , as revealed from 
the numerical simulations. 

After the matter-wave soliton is created, the delocaliz- 
ing transition can be induced by decreasing the strength 
of the OL through the intensity of laser wave, or reduc- 
ing the coefficient of nonlinearity by changing the atomic 
scattering length via the Feshbach resonance. The sig- 
nature of the phenomenon is sudden disintegration of 
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the soliton at the critical point, therefore it will be ac- 
companied by dropping of the atomic density, which can 
be detected by imaging techniques. The resolution of 
absorption-imaging systems currently in use (which is 
~ 7 lira [25]) is sufficient for controlling the properties 
of BEC in the range of a unit cell of the OL. 

Now let us estimate the dimensionless parameters in 
physical units. In typical experiments to date the rele- 
vant parameters arc given by uq = 10^*^ m~"^, — 5.4 
nm, and fcz, = 2n/X = 8.06 • lO'^ m'^ for Rb [27], and 
no = 3-10^^ m-3, = 2.65 nm, and fcz, = 1.07-10'' m"^ 
for Na [28] . The strength of nonlinear atomic interaction 
is given by x = Simoag/k^ = 0.21 for Rb, and x = ^■'^^ 
for Na. Higher or lower values of x may be achieved by 
changing the density of the condensate r?.o, atomic scat- 
tering length fls, or k^. All three parameters can be 
changed independently. The strengths of OLs considered 
above are in the range of usual experimental conditions 
[29] e = 0-=-20 in units of recoil energy Erec = ?i^fc£/(2m). 
Our numerical simulation times t ^ 1000, well satisfying 
the adiabaticity condition, correspond to ~ 0.5 s. The 
average lifetime of BEC in optical trap, limited by three- 
body collisions and off-resonant scattering of lattice pho- 
tons is more than 3 s [30]. Therefore, observation of 
the delocalizing transition of solitons is possible in the 
present experimental conditions. 



V. CONCLUSIONS 

We have studied the existence and delocalizing transi- 
tion of multidimensional solitons in periodic potentials by 
means of the VA and direct numerical integration of the 
Gross-Pitaevskii equation. VA provides the initial con- 
figuration for the soliton parameters which can be used 
in PDE to generate a stable multidimensional soliton. 

The most interesting property of multidimensional (2D 
and 3D) solitons observed in this study is the delocalizing 
transition, which is manifested as irreversible disintegra- 
tion of the soliton at some critical strength of the peri- 
odic potential, or coefficient of nonlincarity. Contrarily, 
ID solitons do not exhibit the delocalizing transition, re- 
taining their integrity over the whole range of parameter 
variations. 

We proposed a physical mechanism for delocalizing 
transition of solitons, according to which the existence 
of a soliton is associated with the existence of quantum 
bound states in the effective potential created by the soli- 
ton. At the point of delocalizing transition, the effective 
potential has the critical strength to support only a sin- 
gle bound state. As the strength becomes weaker than 
the critical value, the effective potential cannot support 
a bound state, and as a consequence, the soliton irre- 
versibly disintegrates. Using the exactly solvable Poschl- 
Tcller potential as approximation for the effective poten- 
tial created by the soliton, we analytically determined 
the existence region for 2D solitons. Numerical simula- 



tions of the Gross-Pitaevskii equation have confirmed the 
validity of the proposed model. 

Although we considered the problem with the emphasis 
on Bose-Einstein condensates in OLs, multidimensional 
solitons and the phenomenon of delocalizing transition 
can be observed in other relevant systems, e.g. optically 
induced nonlinear waveguide arrays [31]. 
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